Robust quantum parameter estimation: coherent magnetometry with feedback 
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We describe the formalism for optimally estimating and controlling both the state of a spin 
ensemble and a scalar magnetic field with information obtained from a continuous quantum limited 
measurement of the spin precession due to the field. The full quantum parameter estimation model 
is reduced to a simplified equivalent representation to which classical estimation and control theory 
is applied. We consider both the tracking of static and fluctuating fields in the transient and steady 
state regimes. By using feedback control, the field estimation can be made robust to uncertainty 
about the total spin number. 
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I. INTRODUCTION 

As experimental methods for manipulating physical 
systems near their fundamental quantum limits improve 
0, 13, 0, 0, , the need for quantum state and param- 
eter estimation methods becomes critical. Integrating 
a modern perspective on quantum measurement theory 
with the extensive methodologies of classical estimation 
and control theory provides new insight into how the lim- 
its imposed by quantum mechanics affect our ability to 
measure and control physical systems 0, S IS • 

In this paper, we illustrate the processes of state es- 
timation and control for a continuously-observed, coher- 
ent spin ensemble (such as an optically pumped cloud 
of atoms) interacting with an external magnetic field. 
In the situation where the magnetic field is either zero 
or well-characterized, continuous measurement (e.g., via 
the dispersive phase shift or Faraday rotation of a far 
off-resonant probe beam) can produce a spin-squeezed 
[Toj state conditioned on the measurement record [TH . 
Spin-squeezing indicates internal entanglement between 
the different particles in the ensemble 12] and promises 
to improve precision measurements 13]. When, however, 
the ambient magnetic environment is either unknown or 
changing in time, the external field can be estimated by 
observing L armor precession in the measurement signal 
0, 0, EE 5 see Figure ^ Recently, we have shown 
that uncertainty in both the magnetic field and the spin 
ensemble can be simultaneously reduced through contin- 
uous measurement and adequate quantum filtering p^ . 

Here, we expand on our recent results 17] involving 
Heisenberg-limited magnetometry by demonstrating the 
advantages of including feedback control in the estima- 
tion process. Feedback is a ubiquitous concept in classi- 
cal applications because it enables precision performance 
despite the presence of potentially large system uncer- 
tainty. Quantum optical experiments are evolving to the 
point where feedback can been used, for example, to sta- 
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bilize atomic motion within optical lattices [J and high 
finesse cavities |5| . Recently, we demonstrated the use of 
feedback on a polarized ensemble of laser-cooled Cesium 
atoms to robustly estimate an applied magnetic field . 
In this work, we investigate the theoretical limits of such 
an approach and demonstrate that an external magnetic 
field can be measured with high precision despite sub- 
stantial ignorance of the size of the spin ensemble. 

The paper is organized as follows. In Section ^ we 
provide a general introduction to quantum parameter 
estimation followed by a specialization to the case of a 
continuously measured spin ensemble in a magnetic field. 
By capitalizing on the Gaussian properties of both coher- 
ent and spin-squeezed states, we formulate the parameter 
estimation problem in such a way that techniques from 
classical estimation theory apply to the quantum system. 
Section Hill presents basic filtering and control theory in a 
pedagogical manner with the simplified spin model as an 
example. This theory is applied in Section HVl where we 
simultaneously derive mutually dependent magnetome- 
try and spin-squeezing limits in the ideal case where the 
observer is certain of the spin number. We consider the 
optimal measurement of both constant and fluctuating 
fields in the transient and steady state regimes. Finally, 
we show in Section that the estimation can be made 
robust to uncertainty about the total spin number by 
using precision feedback control. 



II. QUANTUM PARAMETER ESTIMATION 

First, we present a generic description of quantum 
parameter estimation S S 0- This involves de- 
scribing the quantum system with a density matrix and 
our knowledge of the unknown parameter with a clas- 
sical probability distribution. The objective of parame- 
ter estimation is then to utilize information gained about 
the system through measurement to conditionally update 
both the density matrix and the parameter distribution. 
After framing the general case, our particular example of 
a continuously measured spin ensemble is introduced. 




FIG. 1: (A) A spin ensemble is initially prepared in a coherent state polarized along x, with symmetric variance in the y 
and z directions. Subsequently, a field along y causes the spin to rotate as the z-component is continuously measured. (B) 
Experimental schematic for the measurement process. A far off resonant probe beam traverses the sample and measures the 
z-component of spin via Faraday rotation. The measurement strength could be improved by surrounding the ensemble with 
a cavity. (C) Experimental apparatus subsumed by the Plant block, which serves to map the total field to the photocurrent, 
h{t) ^ y{t). 



A. General problem 

The following outline of the parameter estimation pro- 
cess could be generalized to treat a wide class of prob- 
lems (discrete measurement, multiple parameters), but 
for simplicity, we will consider a continuously measured 
quantum system with scalar Hamiltonian parameter 
and measurement record y{t). 

Suppose first that the observer has full knowledge of 
the parameter 6. The proper description of the system 
would then be a density matrix pe{t) conditioned on the 
measurement record y{t). The first problem is to find a 
rule to update this density matrix with the knowledge 
obtained from the measurement. As in the problem of 
this paper, this mapping may take the form of a stochas- 
tic master equation (SME). The SME is by definition a 
filter that maps the measurement record to an optimal 
estimate of the system state. 

Now if we allow for uncertainty in 0, then a particularly 
intuitive choice for our new description of the system is 

p{t) = I pe{t)p{0,t)dO (1) 
Je 

where p{0^ t) is a probability distribution representing 
our knowledge of the system parameter. In addition to 
the rule for updating each pe{t), we also need to find a 
rule for updating p{6^ t) according to the measurement 
record. By requiring internal consistency, it is possible 
to find a Bayes rule for updating p{0^t) j^. These two 
update rules in principle solve the estimation problem 
completely. 

Because evolving p{t) involves performing calculations 



with the full Hilbert space in question, which is often 
computationally expensive, it is desirable to find a re- 
duced description of the system. Fortunately, it is often 
possible to find a closed set of dynamical equations for 
a small set of moments of p{t). For example, if c is an 
observable, then we can define the estimate moments 

(c)(t) = Tv[p{t)c] 
{Ac'm ^ Tr[p(t)(c-(c))2] 

{0){t) = Jp{e,t)0d0 

{A0^){t) ^ I p{e,t){e-{0)fd0 

and derive their update rules from the full update rules, 
resulting in a set of ?/(t)-dependent differential equations. 
If those differential equations are closed, then this re- 
duced description is adequate for the parameter estima- 
tion task at hand. This situation (with closure and Gaus- 
sian distributions) is to be expected when the system is 
approximately linear. 



B. Continuously measured spin system 

This approach can be applied directly to the problem of 
magnetometry considered in this paper. The problem can 
be summarized by the situation illustrated in Figure ^ 
a spin ensemble of possibly unknown number is initially 
polarized along the x-axis (e.g., via optical pumping), an 
unknown possibly fluctuating scalar magnetic field b di- 
rected along the y-axis causes the spins to then rotate 
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within the x-z plane, and the z-component of the col- 
lective spin is measured continuously. The measurement 
can, for example, be implemented as shown, where we 
observe the difference photocurrent, y{t)^ in a polarime- 
ter which measures the Faraday rotation of a linearly 
polarized far off resonant probe beam travelling along z 
[3,0,^3. The goal is to optimally estimate b{t) via the 
measurement record and unbiased prior information. If 
a control field u{t) is included, as it will be eventually, 
the total field is represented by h{t) = h{t) + u{t). 

In terms of our previous discussion, we have here the 
observable c = \fMJz , where M is the measurement rate 
(defined in terms of probe beam parameters), and the 
parameter = b. When b is known, our state estimate 
evolves by the stochastic master equation pj| 



dpb{t) = -i[H{b), pb{t)]dt ^V[VMJMt)dt 
+V^H[VmJ,] (2^/M^[y{t)dt - {J,)bdt]) pb{t) (2) 

where H{b) = jJyb^ 7 is the gyromagnetic ratio, and 

T>[c]p = cpc^ — {c^cp-\- pc^c)/2 
H[c]p = cp^ pc^ -Tt[{c^ c^)p]p 

The stochastic quantity 2y/Mr][y{t)dt — {Jz)b{t)dt] = 
dW{t) is a Wiener increment (Gaussian white noise with 
variance dt) by the optimality of the filter. The sensi- 
tivity of the photodetection per a/Hz is represented by 
l/2y/Mr], where the quantity r] represents the quantum 
efficiency of the detection. If 7^ = 0, we are essentially ig- 
noring the measurement result and the conditional SME 
becomes a deterministic unconditional master equation. 
If 7^ = 1, the detectors are maximally efficient. Note that 
our initial state p{0) = pb{0) is made equal to a coherent 
state (polarized in x) and is representative of our prior 
information. 

It can be shown that the unnormalized probability 
p(6, t) evolves according to 6] 



dp{b,t) = 4Mr]{J,)b{t)p{b,t)y{t)dt 



(3) 



The evolution Eqs. Q and ^ together with Eq. Q solve 
the problem completely, albeit in a computationally ex- 
pensive way. Clearly, for large ensembles it would be 
advantageous to reduce the problem to a simpler descrip- 
tion. 

If we consider only the estimate moments {Jz)(t)^ 
(AJ^)(t), {b){t), and (A6^)(t) and derive their evolution 
with the above rules, it can be shown that the filtering 
equations for those variables are closed under certain ap- 
proximations. First, the spin number J must be large 
enough that the distributions for Jy and Jz are approx- 
imately Gaussian for an x-polarized coherent state. Sec- 
ond, we only consider times t <C 1/M because the total 
spin becomes damped by the measurement at times com- 
parable to the inverse of the measurement rate. 

Although this approach is rigorous and fail-safe, the re- 
sulting filtering equations for the moments can be arrived 



at in a more direct manner as discussed in Appendix IXI 
Essentially, the full quantum mechanical mapping from 
h(t) to y{t) is equivalent to the mapping derived from 
a model which appears classical, and assumes an actual, 
but random, value for the z component of spin. This 
correspondence generally holds for a stochastic master 
equation corresponding to an arbitrary linear quantum 
mechanical system with continuous measurement of ob- 
servables that are linear combinations of the canonical 
variables 20]. 

From this point on we will only consider the simpli- 
fied Gaussian representation (used in the next section) 
since it allows us to apply established techniques from 
estimation and control theory. The replacement of the 
quantum mechanical model with a classical noise model 
is discussed more fully in the appendix. Throughout 
this treatment, we keep in mind the constraints that 
the original model imposed. As before, we assume J 
is large enough to maintain the Gaussian approximation 
and that time is small compared to the measurement in- 
duced damping rate, t <C 1/M. Also, the description of 
our original problem demands that (AJ^)(0) = J/2 for 
a coherent state [231 • Hence our prior information for 
the initial value of the spin component will always be 
dictated by the structure of Hilbert space. 



III. OPTIMAL ESTIMATION AND CONTROL 

We now describe the dynamics of the simplified rep- 
resentation. Given a linear state-space model (L), a 
quadratic performance criterion (Q) and Gaussian noise 
(G), we show how to apply standard LQG analysis to 
optimize the estimation and control performance 21]. 

The system state we are trying to estimate is repre- 
sented by 



State 



-z{t) 

m 



(4) 



where z{t) represents the small z-component of the col- 
lective angular momentum and b{t) is a scalar field along 
the y axis. 

Our best guess of x(t), as we filter the measurement 
record, will be denoted 



Estimate 



m{t) 



(5) 



As stated in the appendix, we implicitly make the as- 
sociations: z{t) = {Jz){t) = TT[p{t)Jz] and b{t) = 
J p{b^ t)b db^ although no further mention of p{t) or p(6, t) 
will be made. 

We assume the measurement induced damping of J 
to be negligible for short times (Jexp[— Mt/2] ^ J \i 
t <C 1/M) and approximate the dynamics as 
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Dynamics 

dx{t) 

A 
B 



Ax(t)dt + ^u{t)dt 



J J' 
-76 



dWi (6) 



7J 




5]o = 
5^1 = 











(JbF 



where the initial value x{0) for each trial is drawn ran- 
domly from a Gaussian distribution of mean zero and 
covariance matrix XIq. The initial field variance a^o is 
considered to be due to classical uncertainty, whereas the 
initial spin variance a^o is inherently non-zero due to the 
original quantum state description. Specifically, we im- 
pose azo = (AJ^)(0). The Wiener increment dWi{t) has 
a Gaussian distribution with mean zero and variance dt. 
I]i represents the covariance matrix of the last vector in 
Eq. ®. 

We have given ourselves a magnetic field control input, 
u{t)^ along the same axis, y, of the field to be measured, 
b{t). We have allowed b{t) to fiuctuate via a damped 
diffusion (Ornstein-Uhlenbeck) process 



db{t) 



-Ihbit) dt 



The h{t) fiuctuations are represented in this particular 
way because Gaussian noise processes are amenable to 
LQG analysis. The variance of the field at any particu- 
lar time is given by the expectation a^Free = E[6(t)^] = 
cr5i?/276. (Throughout the paper we use the notation 
E[x(t)] to represent the average of the generally stochas- 
tic variable x(t) at the same point in time, over many tra- 
jectories.) The bandwidth of the field is determined by 
the frequency 75 alone. When considering the measure- 
ment of fiuctuating fields, a valid choice of prior might be 
c^50 = cTfoFree, but wc choosc to let (jfoo remain indepen- 
dent. For constant fields, we set (JhFree = 0, but a^o 7^ 0. 

Note that only the small angle limit of the spin motion 
is considered. Otherwise we would have to consider dif- 
ferent components of the spin vector rotating into each 
other. The small angle approximation would be invalid 
if a field caused the spins to rotate excessively, but using 
adequate control ensures this will not happen. Hence, we 
use control for essentially two reasons in this paper: first 
to keep our small angle approximation valid, and, second, 
to make our estimation process robust to our ignorance 
of J. The latter point will be discussed in Section Ivl 

Our measurement of z is described by the process 

Measurement 



y{t)dt Cx(t)dt ^ ./a^dW2{t) 
C ^ [1 0] 

5^2 = CTM = 1/4M7/ 



(8) 



where the measurement shot-noise is represented by the 
Wiener increment dW2{t) of variance dt. Again, ^Jcfm 
represents the sensitivity of the measurement, M is the 
measurement rate (with unspecified physical definition 
in terms of probe parameters), and r] is the quantum 
efficiency of the measurement. The increments dWi and 
dW2 are uncorrelated. 

Following 21], the optimal estimator for mapping y{t) 
to m{t) takes the form 



Estimator 



drn{t) 



m(0) 

Ko(t) 
dJ:{t) 
dt 



Am{t)dt + 'Bu{t)dt 

+Ko(t)[^(t)-Cm(t)]dt 

0' 


^(t)C^^^' 

El + A5](t) + 5](t)A^ 

-5](t)c^x;2"^cx;(t) 

cTcnit) abR{t) 

CTzO 
(750 



(9) 



(10) 
(11) 
(12) 



Eq. (O is the Kalman filter which depends on the 
solution of the matrix Riccati Eq. (|10|) . The Riccati 
equation gives the optimal observation gain Ko(t) for 
the filter. The estimator is designed to minimize the 
average quadratic estimation error for each variable: 
E[{z{t) - z{t))^] and E[{b{t) - b{t))^]. If the model is 
correct, and we assume the observer chooses his prior in- 
formation 5](0) to match the actual variance of the initial 
data Eq, then we have the self-consistent result: 

azE{t) = E[{z{t)-z{t)f]=azR{t) 
a,E{t) = E[{b{t)-b{t)f]=a,R{t) 

Hence, the Riccati equation solution represents both the 
observer gain and the expected performance of an opti- 
mal filter using that same gain. 

Now consider the control problem, which is in many 
respects dual to the estimation problem. We would like 
to design a controller to map y{t) to u{t) in a manner 
that minimizes the quadratic cost function 







Minimized Cost 

H 



P 

Q 



[x^{t)'Px{t) ^u{t)Qu{t)] dt 



(T)Pix(T) 




(13) 



where Pi is the end-point cost. Only the ratio p/q ever 
appears, of course, so we define the parameter A = ^/p/q 
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and use it to represent the cost of control. By setting 
A ^ oo, as we often choose to do in the subsequent anal- 
ysis to simplify results, we are putting no cost on our 
control output. This is unrealistic because, for exam- 
ple, making A arbitrarily large implies that we can apply 
transfer functions with finite gain at arbitrarily high fre- 
quencies, which is not experimentally possible. Despite 
this, we will often consider the limit A ^ oo to set bounds 
on achievable estimation and control performance. The 
optimal controller for minimizing Eq. is 



Total State Dynamics 



Controller 



u{t) 

Kc(t) 
dV{T) 
dT 

v(r = 0) 



-Kc(t)m(t) 
Q-^B^V(r-t) 

P + A^V(T) + V(T)A 

-V(T)BQ-^B^V(T) 
Pi 



(14) 



Here V(T) is solved in reverse time T, which can be in- 
terpreted as the time left to go until the stopping point. 
Thus if T oo, then we only need to use the steady 
state of the V Riccati Eq. (|T5j) to give the steady state 
controller gain Kc for all times. In this case, we can 
ignore the (reverse) initial condition Pi because the con- 
troller is not designed to stop. Henceforth, we will make 
Kc equal to this constant steady state value, such that 
the only time varying coefficients will come from Ko(t). 

In principle, the above results give the entire solution 
to the ideal estimation and control problem. However, 
in the non-ideal case where our knowledge of the system 
is incomplete, e.g. J is unknown, our estimation per- 
formance will suffer. Notation is now introduced which 
produces trivial results in the ideal case, but is helpful 
otherwise. Our goal is to collect the above equations into 
a single structure which can be used to solve the non- 
ideal problem. We define the total state of the system 
and estimator as 



Total State 



e{t) 



'xit) 
m{t) 



■z{t) 

m 
m 
m 



(16) 



Consider the general case where the observer assumes 
the plant contains spin J', which may or may not be equal 
to the actual J. All design elements depending on J' in- 
stead of J are now labelled with a prime. Then it can 
be shown that the total state dynamics from the above 
estimator-controller architecture are a time-dependent 
Ornstein-Uhlenbeck process 



de{t) 

a{t) 

m 



a{t)0{t)dt^ (3{t)dW{t) 
A -BK[ 



(17) 



c 



■0 
















_0 






BKc 

0' 




where the covariance matrix of dW is dt times the iden- 
tity. Now the quantity of interest is the following covari- 
ance matrix 



Total State Covariance 



0(t) ^ mt)o^{t)\ 







C^zb 


C^zz 


^zh 


(15) 


C^zb 


C^bb 


CTbz 


^bb 




Cfzz 


Cfbz 


Cfzz 


^zb 




-^zb 


^bb 


^zb 


^66^ 



(18) 



(^ZZ 

O'zb 



E[z{tf] 
E[zmt)] 



It can be shown that this total covariance matrix obeys 
the deterministic equations of motion 



Total State Covariance Dynamics 

dS{t) 



dt 



a{t)e{t)^e{t)a' {t) 
rt 

00 exp 



exp 



- / a{t')dt' 
Jo 

dt' exp 



a{s)ds 



X exp 



{s)ds 



mp^it) (19) 

- f a^{t')dt' 
. Jo 

(20) 



©0 



'(Jzo 0- 

(760 



0. 



Eq. H2()|l is the matrix form of the standard integrating 
factor solution for time-dependent scalar ordinary differ- 
ential equations |23| • Whether we solve this problem nu- 
merically or analytically, the solution provides the quan- 
tity that we ultimately care about 

Average Magnetometry Error 

(^bE{t) 



E[m-h{t)f\ 

E[h\t)]+E[h\t)]-2-E[h{t)h{t)] 



(21) 
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When all parameters are known (and J' = J), this 
total state description is unnecessary because (JbE{t) = 
(Ti)R{t). This equality is by design. However, when the 
wrong parameters are assumed (e.g., J' ^ J) the equality 
does not hold (JbE{t) 7^ crbnit) and either Eq. (|T9j) or Eq. 
(^nj) must be used to find (JhEit)- Before addressing this 
problem, we consider in detail the performance in the 
ideal case, where all system parameters are known by 
the observer, including J. 

At this point, we have defined several variables. For 
clarity, let us review the meaning of several before contin- 
uing. Inputs to the problem include the field fiuctuation 
strength ct^f, Eq. 0, and the measurement sensitivity 
ctm, Eq. 0. The prior information for the field is la- 
belled cr^o, Eq. (|T^ . The solution to the Riccati equa- 
tion is (Thnif)^ Eq. (|TH) . and is equal to the estimation 
variance (JbE{t)^ Eq. (|21|) . when the estimator model is 
correct. In the next section, we additionally use (755, Eq. 

and (JhT{t)i Eq. ((^ . to represent the steady state 
and transient values of (JbE{t) respectively. 



IV. OPTIMAL PERFORMANCE: J KNOWN 

We start by observing qualitative characteristics of the 
6-estimation dynamics. Figure [21 shows the average esti- 
mation performance, (TbR{t)^ as a function of time for a 
realistic set of parameters. Notice that a^R is constant 
for small and large times, below ti and above ^2- If (Jbo is 
non-infinite then the curve is constant for small times, as 
it takes some time to begin improving the estimate from 
the prior. If a^o is infinite, then ti = and the sloped 
transient portion extends towards infinity as t ^ 0. At 
long times, a^R will become constant again, but only if 
the field is fiuctuating (ct^f 7^ and 75 7^ 0). The per- 
formance saturates because one can track a field only so 
well if the field is changing and the signal-to-noise ratio is 
finite. If the field to be tracked is constant, then ^2 = 00 
and the sloped portion of the curve extends to zero as 
t ^ 00 (given the approximations discussed in Section 
IIIB|) . After the point where the performance saturates 
{t ^ ^2), all of the observer and control gains have be- 
come time-independent and the filter can be described 
by a transfer function. 

However, as will be shown, applying only this steady 
state transfer function is non-optimal in the transient 
regime {ti <C t <C ^2), because the time dependence of the 
gains is clearly crucial for optimal transient performance. 



A. Steady state performance 

We start by examining the steady state performance 
of the filter. At large enough times (where we have yet to 
define large enough), Kq becomes constant and if we set 
T ^ 00 (ignoring the end point cost), then Kc is always 
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FIG. 2: The Riccati equation solution gives the ideal field 
estimation performance. The parameters used here are J — 
10^ azo = J/2 (for ensemble of spin-l/2's), 7 = 10^ M = 
10^, cTfoo = cTbFree = 1- The solutiou staits at the free field 
fluctuation variance and saturates at ats- The plot is not 
valid at times t > 1/M. 



constant. Setting dT,/dt = and d\/dt = we find: 




■76) 



[A 1/(1 



Ko(t) 

Kc(t) 
where A = 

Now assuming the gains to be constant, we can derive 
the three relevant transfer functions from y{t) to rn{t) 
{z and b) and u. We proceed as follows. First, we ex- 
press the estimates in terms of only themselves and the 
photocurrent 



dm{t) 
dt 



Am(t) +Bii(t) +Ko(^(t) - Cm(t)) 

= Am{t) + B(-Kcm(t)) + Ko(^(t) - Cm{t)) 
= {A-BKc-KoC)m{t)^Koy{t) 



To get the transfer functions, we take the Laplace trans- 
form of the entire equation, use differential transform 
rules to give s factors (where s = juo, j = a/— T), ignore 
initial condition factors, and rearrange terms. However, 
this process only gives meaningful transfer functions if 
the coefficients Kq and Kc are constant. Following this 
procedure, we have 



m{s) 



u{s) 



(5l - A + BKc + KoC)-^Koy{s) 

Gm{s)yis) 



= -Kcm{s) 
= -Kc{s-A- 
= Gu{s)y{s) 



BK 



c 



-KoC)-'Koy{s) 
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where 



Gm{s) 



The three transfer functions {Gz{s)^ Gb{s), and Gu{s)) 
serve three different tasks. If estimation is the concern, 
then Gb{s) wih perform optimahy in steady state. No- 
tice that, while the Riccati solution is the same with and 
without control (Kc non-zero or zero), this transfer func- 
tion is not the same in the two cases. So, even though 
the transfer functions are different, they give the same 
steady state performance. 

Let us now consider the controller transfer function 
Gu{s) in more detail. We find the controller to be of the 
form 



Gu{s) = Gu,DC 



1 + s/iJH 



1 + (1 + s/uq)s/uol 



(22) 



Here each frequency uo represents a transition in the Bode 
plot of Figure |2l A similar controller transfer function is 
derived via a different method in Appendix O 

If we are not constrained experimentally, we can make 

the approximations ^ \J (JhF I I'^lJ and ^ 

ll\/(^M/crhF giving 



Gu{s) 

LUC 
Gu.DC 



Gu,DC 



s/ujh 



1 + s/ljl 




= 2u;h 



Gu 



Gu.DC—^ = —\l — 




where Gn,AC is the gain at high frequencies {uj > uoh) 
and we find the closing frequency uoc from the condition 
\Pz{j^c)Gu{j^c)\ = 1, with the plant transfer function 
being the normal integrator Pz{s) = jJ/s. Notice that 
the controller closes in the very beginning of the flat high 
frequency region (hence with adequate phase margin) be- 
cause ujc = 2u;H' 

Finally, consider the steady state estimation perfor- 
mance. These are the same with and without control 
(hence A-independent) and, under the simplifying as- 
sumption 7 J 75 y^cTM/crbF^ are given by 
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FIG. 3: The Bode plot of Gu{s), the transfer function of the 
filter in steady state, for a typical parameter regime. Notice 
that the controller closes the plant with adequate phase mar- 
gin to avoid closed-loop instability. At high frequencies the 
controller rolls off at uuq if A ^ oo. 



If the estimator reaches steady state at t <C then 
the above variance azR represents a limit to the amount 
of spin-squeezing possible in the presence of fluctuating 
fields. 

Also the J scaling of the saturated field sensitivity 
o'bR J~^/^ is not nearly as strong as the J scaling in 
the transient period (JbR oc J~^. Next, we demonstrate 
this latter result as we move from the steady state analy- 
sis to calculating the estimation performance during the 
transient period. 



B. Transient performance 

We now consider the transient performance of the ideal 
filter: how quickly and how well the estimator-controller 
will lock onto the signal and achieve steady state per- 
formance. In many control applications, the transient 
response is not of interest because the time it takes to 
acquire the lock is negligible compared to the long steady 
state period of the system. However, in systems where 
the measurement induces continuous decay, this transient 
period can be a significant portion of the total lifetime of 
the experiment. 

We will evaluate the transient performance of two dif- 
ferent filters. First, we look at the ideal dynamic ver- 
sion, with time dependent observer gains derived from 
the Riccati equation. This limits to a transfer function 
at long times when the gains have become constant. Sec- 
ond, we numerically look at the case where the same 
steady state transfer functions are used for the entire 
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TABLE I: Field tracking error, abR{t), for different initial variances of b and z. 
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duration of the measurement. Because the gains are not 
adjusted smoothly, the small time performance of this es- 
timator suffers. Of course, for long times the estimators 
are equivalent. 



1. Dynamic estimation and control 

Now consider the transient response of (giving 
Ko(t)). We will continue to impose that V (thus Kc) is 
constant because we are not interested in any particular 
stopping time. 

The Riccati equation for 5](t) (Eq. COJ) appears dif- 
ficult to solve because it is non-linear. Fortunately, it 
can be reduced to a much simpler linear problem. See 
Appendix IbI for an outline of this method. 

The solution to the fluctuating field problem (ct^f / 
and 76 7^ 0) is represented in Figure El This solution is 
simply the constant field solution {cFhF = and 75 = 0) 
smoothly saturating at the steady state value of Eq. ((211) 
at time ^2- Thus, considering the long time behavior of 
the constant field solution will tell us about the transient 
behavior when measuring fluctuating fields. Because the 
analytic form for the constant field solution is simple, we 
consider only it and disregard the full analytic form of 
the fluctuating field solution. 

The analytic form of 5](t) is highly instructive. The 
general solutions to (JbR{t) and (JzR{t)^ with arbitrary 
prior information a^o and a^o, are presented in the cen- 
tral entries of Tables ^ and ^ respectively. The other 
entries of the tables represent the limits of these some- 
what complicated expressions as the prior information 
assumes extremely large or small values. Here, we notice 
several interesting trade-oflFs. 

First, the left hand column of Tabled is zero because 
if a constant field is being measured, and we start with 



complete knowledge of the field (0-50 = 0), then our job 
is completed trivially. Now notice that if a^o and a^o are 
both non-zero, then at long times we have the lower right 
entry of Table U 



CTbRit) 



cTbrit) 



(25) 



This is the same result one gets when the estimation pro- 
cedure is simply to perform a least-squares line fit to the 
noisy measurement curve for constant fields. (Note that 
all of these results are equivalent to the solutions of ^3 5 
but without J damping.) If it were physically possible 
to ensure a^o = 0, then our estimation would improve 
by a factor of four to the upper right result. However, 
quantum mechanics imposes that this initial variance is 
non-zero (e.g., azo = J/2 for a coherent state and less, 
but still non-zero, for a squeezed state), and the upper 
right solution is unattainable. 

Now consider the dual problem of spin estimation per- 
formance (JzR(t) as represented in Table Ull where we can 
make analogous trade-off observations. If there is no field 
present, we set cr^o = and 



CTzRif) 



O'zOO'M 



Cm -\- ta 



zO 



(26) 



When cFzRit) is interpreted as the quantum variance 
(AJ^)(t), this is the ideal (non-damped) conditional spin- 
squeezing result which is valid at t <C 1/M, before damp- 
ing in J begins to take effect 19]. If we consider the 
solution for t ^ 1/JM, we have the lower left entry of 
Table ITll (7zR{t) = crM/t- However, if we must include 
constant field uncertainty in our estimation, then our es- 
timate becomes the lower right entry (TzR{t) = 4.(7 m/^ 
which is, again, a factor of four worse. 

If our task is field estimation, intrinsic quantum me- 
chanical uncertainty in z limits our performance just as. 
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if our task is spin-squeezed state preparation, field un- 
certainty limits our performance. 



2. Transfer function estimation and control 

Suppose that the controller did not have the capability 
to adjust the gains in time as it tracked a fluctuating field. 
One approach would then be to apply the steady state 
transfer functions derived above for the entire measure- 
ment. While this approach performs optimally in steady 
state, it approaches the steady state in a non-optimal 
manner compared to the dynamic controller. Figure ^ 
demonstrates this poor transient performance for track- 
ing fluctuating fields of differing bandwidth. Notice that 
the performance only begins to improve around the time 
that the dynamic controller saturates. 

Also notice that the transfer function Gb{s) is depen- 
dent on whether or not the state is being controlled, i.e. 
whether or not A is zero. The performance shown in Fig- 
ure ^ is for one particular value of A, but others will give 
different estimation performances for short times. Still, 
all of the transfer functions generated from any value of 
A will limit to the same performance at long times. Also, 
all of them will perform poorly compared to the dynamic 
approach during the transient time. 



V. ROBUST PERFORMANCE: J UNKNOWN 

Until this point, we have assumed the observer has 
complete knowledge of the system parameters, in particu- 
lar the spin number J. We will now relax this assumption 
and consider the possibility that, for each measurement, 
the collective spin J is drawn randomly from a particu- 
lar distribution. Although we will be ignorant of a given 
J, we may still possess knowledge about the distribution 
from which it is derived. For example, we may be certain 
that J never assumes a value below a minimal value Jmin 
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or above a maximal value 



This is a realistic ex- 



perimental situation, as it is unusual to have particularly 
long tails on, for example, trapped atom number distri- 
butions. We do not explicitly consider the problem of J 
fluctuating during an individual measurement, although 
the subsequent analysis can clearly be extended to this 
problem. 

Given a J distribution, one might imagine completely 
re-optimizing the estimator-controller with the full dis- 
tribution information in mind. Our initial approach is 
more basic and in line with robust control theory: we 
design our filter as before, assuming a particular J', then 
analyze how well this filter performs on an ensemble with 
J 7^ J^ With this information in mind, we can decide if 
estimator-controllers built with J' are robust, with and 
without control, given the bounds on J. We will find 
that, under certain conditions, using control makes our 
estimates robust to uncertainty about the total spin num- 
ber. 
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FIG. 4: Estimation performance for estimators based on 
the dynamic gain solution of the Riccati equation, compared 
against estimators with constant estimation gain. The latter 
are the transfer function limits of the former, hence they have 
the same long term performance. Three different bandwidth 
b processes are considered. 



The essential reason for this robustness is that when 
a control field is applied to zero the measured signal, 
that control field must be approximately equal to the 
field to be tracked. Because J is basically an effective 
gain, variations in J will affect the performance, but not 
critically, so the error signal will still be approximately 
zero. If the applied signal is set to be the estimate, then 
the tracking error must also be approximately zero. (See 
Appendix [0 for a robustness analysis along these lines 
in frequency space. These methods were used in but 
neglect the transient behavior.) 

Of course, this analysis assumes that we can apply 
fields with the same precision that we measure them. 
While the precision with which we can apply a field is 
experimentally limited, we here consider the ideal case of 
infinite precision. In this admittedly idealized problem, 
our estimation is limited by only the measurement noise 
and our knowledge of J. 

First, to motivate this problem, we describe how poorly 
our estimator performs given ignorance about J without 
control. 



A. Uncontrolled ignorance 

Let us consider the performance of our estimation pro- 
cedure at estimating constant fields when J' ^ J . In 
general, this involves solving the complicated total co- 
variance matrix Eq. (|20|) . However, in the long time limit 
{t ^ 1/JM) of estimating constant fields, the procedure 
amounts to simply fitting a line to the noisy measurement 
with a least-squares estimate. Suppose we record an 
open-loop measurement which appears as a noisy sloped 
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line for small angles of rotation due to the Larmor pre- 
cession. Regardless of whether or not we know J, we can 
measure the slope of that line and estimate it to be m. If 
we knew J, we would know how to extract the field from 
the slope correctly: b = m/jJ. If we assumed the wrong 
spin number, J' ^ J, we would get the non-optimal esti- 
mate: h' = m/^J' = hJ/J'. 

First assume that this is a systematic error and J is 
unknown, but the same, on every trial. We assume that 
the constant field is drawn randomly from the a^o distri- 
bution for every trial. In this case, if we are wrong, then 
we are always wrong by the same factor. It can be shown 
that the error always saturates 

(JhE ^ (1 - /)^Cr5o 

where f = J / J' . Of course, because this error is sys- 
tematic, the variance of the estimate does not saturate, 
only the error. This problem is analogous to ignorance 
of the constant electronic gains in the measurement and 
can also be calibrated away. 

However, a significant problem arises when, on every 
trial, a constant h is drawn at random and J is drawn 
at random from a distribution, so the error is no longer 
systematic. In this case, we would not know whether to 
attribute the size of the measured slope to the size of J or 
to the size of h. Given the same h every trial, all possible 
measurement curves fan out over some angle due to the 
variation in J. After measuring the slope of an individual 
line to beyond this fan-out precision, it makes no sense 
to continue measuring. 

We should also point out procedures for estimating 
fields in open-loop configuration, but without the small 
angle approximation. For constant large fields, we could 
observe many cycles before the spin damped significantly. 
By fitting the amplitude and frequency independently, 
or computing the Fourier transform, we could estimate 
the field somewhat independently of J, which only de- 
termines the amplitude. However, the point here is that 
h might not be large enough to give many cycles before 
the damping time or any other desired stopping time. In 
this case, we could not independently fit the amplitude 
and frequency because they appear as a product in the 
initial slope. Similar considerations apply for the case of 
fiuctuating h and fiuctuating J. See |23|, for a complete 
analysis of Bayesian spectrum analysis with free induc- 
tion decay examples. 

Fortunately, using precise control can make the estima- 
tion process relatively robust to such spin number fiuc- 
tuations. 



B. Controlled ignorance: steady state performance 

We first analyze how the estimator designed with J' 
performs on a plant with J at tracking fiuctuating fields 
with and without control. To determine this we calculate 
the steady state of Eq. ([TUI) . 




Time 



FIG. 5: Steady state estimation performance for estima- 
tor designed with J' — 10^, and actual spin numbers: 
J = f X [0.5, 0.75, 1, 1.25, 2, 10, 100]. Other parameters: 
7 = 10^, M = 10^, 7b = 10^ abFree = 1 (fluctuating field), 
A = 0.1 (this is large enough to satisfy large- A limits discussed 
in text). The inset compares the normalized robust estima- 
tion performance (curve) at a particular time, to the ideal 
performance (line) when J is known. 



For the case of no control (A = 0), we simplify the re- 
sulting expression by taking the same large J' approxima- 
tion as before. This gives the steady state uncontrolled 
error 



CTbE 



(1 



(1 - f)^(TbFr 



where f = J / J' . Because the variance of the fiuctuating 
h is (JbFree^ the Uncontrolled estimation performs worse 
than no estimation at all if /> 2. 

On the other hand, when we use precise control the 
performance improves dramatically. We again simplify 
the steady state solution with the large J' and A assump- 
tions from before, giving 



1 + / 

2/ 
1 + / 
2/ 



3/4 1/4 



CJbs{J') 



where (Jhs{Ji J') is the steady state controlled error when 
a plant with J is controlled with a J' controller and 
^hs{J') is the error when J = J' . One simple interpreta- 
tion of this result is that if we set J' to be the minimum 
of the J distribution (/ > 1) then we never do worse 
than cFi)s{J') and we never do better than twice as well 
(/ oo). See Figure [51 for a demonstration of this per- 
formance. 
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on J' = Jmin and we would be guaranteed at least 
the performance (JbT{Jmin) and at best the performance 

CrbT{Jmin)/^' 

This approach would be suitable for experimental situ- 
ations because typical J distributions are narrow: the dif- 
ference between Jmin and Jmax is rarely greater than an 
order of magnitude. Thus, the overall sacrifice in perfor- 
mance between the ideal case and the robust case would 
be small. The estimation performance still suffers be- 
cause of our ignorance of J, but not nearly as much as in 
the uncontrolled case. 



VI. CONCLUSION 



FIG. 6: Transient estimation performance for controller de- 
signed with J' = 10^, and actual spin numbers: J — J' x 
[0.75, 1, 1.25, 2, 10, 100, 1000]. Other parameters: 7 = 10^ 
M = 10^, 76 = 0, (JbFree = (coustaut field), A = 1. Note 
that this behavior is vaUd for t < 1/M = 10~^. The in- 
set compares the normalized robust estimation performance 
(curve) at a particular time, to the ideal performance (line) 
when J is known. 

C. Controlled ignorance: transient performance 

Now consider measuring constant fields with the wrong 
assumed J^ Again, when control is not used, the error 
saturates at 

When control is used, the transient performance again 
improves under certain conditions. The long time tran- 
sient solution of Eq. (|T^ is difficult to manage analyti- 
cally, yet the behavior under certain limits is again sim- 
ple. For large A and J', and for /> 1/2, we numerically 
find the transient performance to be approximately 

= (^)--(^') (27) 

where cFbriJ^ J') is the transient controlled error when a 
plant with J is controlled with a J' controller and (TbT{J') 
is the error when J = J' . See Figure [SI for a demon- 
stration of this performance for realistic parameters. As 
/ ^ oo the /-dependent pre-factor saturates at a value of 
1/4. However, as / 1/2 then the system takes longer 
to reach such a simple asymptotic form, and the solution 
of Eq. (|27|) becomes invalid. 

Accordingly, one robust strategy would be the follow- 
ing. Suppose that the lower bound of the J-distribution 
was known and equal to Jmin- Also assume that 
^hT^Jmin) represents an acceptable level of performance. 
In this case, we could simply design our estimator based 



The analysis of this paper contained several key steps 
which should be emphasized. Our first goal was to outline 
the proper approach to quantum parameter estimation. 
The second was to demonstrate that reduced representa- 
tions of the full filtering problem are relevant and con- 
venient because, if a simple representation can be found, 
then existing classical estimation and control methods 
can be readily applied. The characteristic that led to 
this simple description was the approximately Gaussian 
nature of the problem. Next, we attempted to present 
basic classical filtering and control methodology in a self- 
contained, pedagogical format. The results emphasized 
the inherent trade-offs in simultaneous estimation of dis- 
tinct, but dynamically coupled, system parameters. Be- 
cause these methods are potentially critical in any field 
involving optimal estimation, we consider the full expo- 
sition of this elementary example to be a useful resource 
for future analogous work. 

We have also demonstrated the general principle that 
precision feedback control can make estimation robust 
to the uncertainty of system parameters. Despite the 
need to assume that the controller produced a precise 
cancellation field, this approach deserves further investi- 
gation because of its inherent ability to precisely track 
broadband field signals 0- It is anticipated that these 
techniques will become more pervasive in the experimen- 
tal community as quantum systems are refined to levels 
approaching their fundamental limits of performance. 
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APPENDIX A: SIMPLIFIED REPRESENTATION 
OF THE PLANT 

In Section ^ we outlined a general approach to quan- 
tum parameter estimation based on the stochastic mas- 
ter equation (SME), but subsequently we derived optimal 
observer and controller gains from an explicit representa- 
tion of the plant dynamics (Eq. 0). This representation 
appears classical in that the plant state is given by a 
scalar variable, rather than a density operator. In this 
Section we present a derivation of this simplified repre- 
sentation and discuss the equivalence of our approach to 
the original quantum estimation problem. 

From the perspective of quantum filtering theory we 
will simply show that a Gaussian approximation to the 
relevant SME can be viewed as a Kalman filter, which in 
turn induces a simplified representation of the dynamics 
of the spin state. In this simplified representation the 
quantum state of the spin system is replaced by a scalar 
variable, and {Jz){t) is viewed as the optimal estimate 
of the random process z{t). Equations for dz{t) and its 
relation to the observed photocurrent y{t)dt are given in 
Eqs. (|A3|) and (|A6|) . which have the convenient property 
of being formally time-invariant. The technical approach 
in the main body of the text is then to replace Eq. (|A1|) . 
which is derived from the SME, by a state-space observer 
derived directly from the simplified model of Eq. (|A3|) . 
By doing so we achieve transparent correspondence with 
classical estimation and control theory. We should note 
that the diagrams in Figure[7|indicate signal fiows and de- 
pendencies in a way that is quite at odds with the quan- 
tum filtering perspective. This Figure is meant solely 
to motivate the simplified model (Eq. ()A3|) ) for readers 
who prefer a more traditional quantum optics perspec- 
tive, in which the Ito increment in the SME corresponds 
to optical shot-noise (as opposed to an innovation pro- 
cess derived from the photocurrent) and the SME itself 
plays the role of a 'physical' evolution equation mapping 
h{t) to y{t)dt. 

Adopting the latter perspective, let us briefiy discuss 
(with reference to the top diagram in Figure^ the overall 
structure of our estimation problem. The physical sys- 
tem that exists in the laboratory (the spins and optical 
probe beam) acts as a transducer, whose key role in the 
magnetometry scheme is to imprint a statistical signature 
of the magnetic field h{t) onto the observable photocur- 
rent, y{t)dt. Hence whatever theoretical model we adopt 
for describing the spin and probe dynamics must pro- 
vide an accurate description of the mapping from h{t) 
to y{t)dt, as represented by the Plant in Figure Qp. An 
open- loop estimator, designed on the basis of this plant 
model, would construct a conditional probability distri- 
bution for h{t) based on passive observation of y{t)dt. In 
a closed-loop estimation procedure we would allow the 
controller to apply compensation fields to the system in 
order to gain accuracy and/or robustness. In either case, 
the essential role of the spin-probe (plant) model in the 
design process is to provide an accurate description of 




FIG. 7: Equivalent models for the filtering problem (see dis- 
cussion at the beginning of Appendix 0. Each version can 
be inserted into the Plant block of Figure QC. The filters all 
presume complete knowledge of h{t) = b{t) + u{t). 



the infiuence of an arbitrary time-dependent field h{t) 
on the photocurrent y{t)dt. Note that the consideration 
of arbitrary h{t) subsumes all possible effects of real-time 
feedback. 

Thomsen and co-workers 19j have derived an accurate 
plant model for our magnetometry problem, in the form 
of an SME (Eq. ((2|)). Following a common convention 
in quantum optics, let us here write this SME and the 
corresponding photocurrent equation in the form 



dp{t) 



y{t)dt 



-i dt[H{h), pit)] + V[VMJ^]p{t) dt 
{J,){t)dt^^/^dW{t), 



where H{h) = jhJy and p{t) is the state of the spin sys- 
tem conditioned on the measurement record y(t)dt. The 
quantity dW{t) is a Wiener increment that heuristically 
represents shot-noise in the photodetection process 24], 
and these are to be interpreted as Ito stochastic differen- 
tial equations. If h{t) and y{t)dt are considered as input 
and output signals, respectively, this pair of equations 
jointly implement a plant transfer function as depicted 
in Figure with p{t) taking on the role of the plant 
state. 

For a large spin ensemble, however, p{t) will have very 
high dimension and it would be impractical to utilize the 
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full SME for design purposes. It is straightforward to 
derive a reduced model by employing a moment expan- 
sion for the observable of interest. Extracting the condi- 
tional expectation values of the first two moments of 
from SME gives the following scalar stochastic differen- 
tial equations: 



dW{t) 



d{J,)(t) = j{J,){t)h{t)dt^ {^J^) 

-ij{[Ajljy]){t)h{t)dt 



If the spins are initially fully polarized along x and 
the spin angle ~ {Jz)/{Jx) is kept small (e.^., by active 
control), then, by using the evolution equation for the x 
component, we can show {Jx){t) ~ Jexp[— Mt/2] ^ J 
for times t < 1/M. Making the Gaussian approxima- 
tion at small times, the third order terms (AJ|) and 
—ij{[AJ'^^Jy]){t)h{t) can be neglected. The Holstein- 
Primakoff transformation [23|, commonly used in the 
condensed matter physics literature, makes it possible 
to derive this Gaussian approximation as an expansion 
in 1/J. Both of the removed terms can be shown to 
be approximately 1 / Ja/J smaller than the retained non- 
linear term. Additionally, the second removed term will 
be reduced if h{t) ~ by active control. 

These approximations give 



d{Jz)(t) = -fJh{t)dt^ 
{AJl)\t) 



ctm 



dt 



dW{t) (Al) 



(A2) 



which constitute a Gaussian, small-time approximation 
to the full SME that represents the essential dynamics 
for magnetometry. Note that we can analytically solve 



(AJ|)(t) 



(AJ|)(0)cTM 

(TM + (Aj2)(0)t' 



where (AJ|)(0) = J/2 for an initial coherent spin state. 

At this point we may note that Eqs. (|A1|) and (|A2|) 
have the algebraic form of a Kalman filter. (This is not 
at all surprising since the SME, as written in Eq. ((SJ, 
represents an optimal nonlinear filter for the reduced spin 
state 9j and our subsequent approximations have en- 
forced both linearity and sufficiency of second-order mo- 
ments.) Viewed as such, the quantity {Jz)(f) would rep- 
resent an optimal (least square) estimate of some under- 
lying variable z{t) based on observation of a signal (i^(t), 
and {AJ^){t) would represent the uncertainty (variance) 
of this estimate. It thus stands to reason that we might 
be able to simplify our magnetometry model even further 
if we could find an 'underlying' model for the evolution 



of z{t) and (i^(t), for which our equations derived from 
the SME would be the Kalman filter. 

It is not difficult to do so, and indeed a very simple 
model suffices: 



dz{t) = -fJh{t)dt 

d^{t) = z{t)dt^ ./o^dW(t) 



(A3) 



where dW(t) is an Wiener increment that is distinct from 
(though related to) dW{t). In order to match initial con- 
ditions with the equations derived from the SME, we 
should assume that the expected value oi z{t = 0) is 
zero and that the variance of our prior distribution for 
z{0) is J/2. Written in canonical form, the Kalman filter 
for this hypothetical system is then 



dzit) 

dcFzRit) 



-fJh{t)dt ■ 



^zR{t) [di{t) - Z{t)dt] 



dt. 



Here z{t) is the optimal estimate of z{t) and (JzR(t) is the 
variance. We exactly recover the SME model, Eqs. (EH) 
and (|A2|) . by the identifications 



[dm 



m 

z{t)dt] 



dW{t). 



(A4) 



It is important to note that the quantity 
[d^{t) — zdt] I ^gm represents the so-called innova- 
tion process of this Kalman filter, and it is thus 
guaranteed (by least-squares optimality of the filter 26]) 
to have Gaussian white-noise statistics. Hence we have 
solid grounds for identifying it with the Ito increment 
appearing in the SME. 

Given this insight, we see that our original magnetom- 
etry problem can equivalently be viewed in a way that 
corresponds to the middle diagram of Figure [7\ In this 
version, we posit the existence of a hidden transducer 
that imprints statistical information about the magnetic 
field h{t) onto a signal d^{t). A Kalman filter receives 
this signal, and from it computes an estimate z{t) as well 
as an innovation process dW{t). (Note that as indicated 
in the diagram, the Kalman filter will only function cor- 
rectly if it 'has knowledge of the true magnetic field h{t) 
in the way that a physical system would, but this is not 
an important point for what follows.) According to the 
model equations, the Kalman filter then emits the fol- 
lowing signal to be received by our photo-detector: 



y{t)dt z{t)dt^ ./^dW{t). 



(A5) 



Note that dW{t) now appears as an internal variable to 
the Kalman filter, computed from the input signal d^^it) 
and the recursive estimate i(t), while the inherent ran- 
domness is referred back to dW(t). Although this may 
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seem like an unnecessarily complicated story, it should 
be noted that the compound model with z{t) and the 
Kalman filter predicts an identical transfer function from 
h{t) to the experimentally-observed signal y{t)dt to that 
of the equations originally derived from the SME (top 
diagram in Figure Ej). Hence, for the purposes of ana- 
lyzing and designing magnetometry schemes, these are 
equivalent models. 

Combining several definitions above we find 



directly from the stochastic master Eq. 



y{t)dt = z(t)dt- 



.[d(,(t) - z(t)dt] 



(A6) 



It thus follows that in the compound model, the Kalman 
filter actually implements a trivial transfer function and 
can in fact be eliminated from the diagram. Doing this, 
we obtain the simplified representation in the bottom dia- 
gram of Figure [3 Here the perspective is to pretend that 
the internal dynamics of the transducing physical system 
corresponds to the simplified model (Eq. (|A3|) ). since we 
can do so without making any error in our description 
of the effect of h{t) on the recorded signal. We thus 
conclude that for the purposes of open- or closed-loop 
estimation of h(t)^ filters and controllers can in fact be 
designed — without loss of performance — using the sim- 
plified model (Eq. (IX3|) ). 

It is interesting to note that z{t) can loosely be in- 
terpreted as a 'classical value' of the spin projection J^. 
Since the operator Jz is a back-action evading observable, 
the continuous measurement we consider is quantum non- 
demolition and its backaction on the system state is 
minimal (conditioning without disturbance). Hence if 
h{t) = 0, we may think of the measurement process as 
gradually 'collapsing' the quantum state of the spin sys- 
tem from an initial coherent state towards an eigenstate 
of Jz] the hidden variable z{t) in the simplified model 
Eq. (|A3|) would then represent the eigenvalue correspond- 
ing to the ultimate eigenstate, and z{t) = {Jz){t) in the 
Kalman filter would be our converging estimate of it. 
(Again, this is as expected from the abstract perspec- 
tive of quantum filtering theory for open quantum sys- 
tems.) Conditional spin-squeezing in this case can then 
be understood as nothing more than the reduction of 
our uncertainty as to the underlying value of z — as 
we acquire information about z through observation of 
d£,{t) = y{t)dt^ our uncertainty (JzR(f) ^ {AJ'^){t) nat- 
urally decreases below its initial coherent-state value of 
J/2. Still, the quantum- mechanical nature of the spin 
system is not without consequence, as it is known that 
continuous QND measurement produces entanglement 
among the spins in the ensemble jl^ . 

It seems worth commenting on the fact that Eq. (|A3|) 
clearly predicts stationary statistics for the photocur- 
rent y{t)dt^ whereas Eq. (|A1|) contains a time-dependent 
diffusion coefficient that might color the statistics of 
y{t)dt = {Jz){t)dt + y/cTM dW{t). In fact there is no 
discrepancy. It is possible 24] to derive the second or- 
der time correlation function of the observed signal y{t)dt 



{y{t)y{t^r)) 



{{Jz{t)Jz{t- 



r))^{Jz{t^r)Jzm/^ 



(This result could also be obtained from the standard 
input-output theory of quantum optics.) Since the mas- 
ter equation results in linear equations for the mean val- 
ues {Jx{t)) and {Jz{t)) the quantum regression theorem 
pT^ allows the correlation functions {J zif) J zit ^ r)) and 
{Jz{t + r)Jz{t)) to be calculated explicitly. In this pa- 
per we are most interested in the early time evolution for 
which we obtain the expressions 



{y{t)) = (Mt)) 

(y(0)y(t))-(2/(0))(2/(t)) 



jbJt 
1 

47]M 
+0{t''). 



5{t) + {Aj!m 



These correlation functions correspond to a white noise 
signal which is a linear ramp with gradient jbJ with a 
random offset of variance (AJ^)(0), in perfect agreement 
with our simplified model Eq. (| A3|) . If the statistics of 
y{t) were Gaussian these first and second order moments 
would be enough to characterize the signal completely, 
and indeed for sufficiently large J the problem does be- 
come effectively Gaussian. 

As a final comment we note that the essential step in 
the above discussion is to observe that the equations for 
the first and second order moments of the quantum state 
derived from the stochastic master equation correspond 
to a Kalman filter for some classical model of a noisy 
measurement. This correspondence holds for the stochas- 
tic master equation corresponding to an arbitrary linear 
quantum mechanical systems with continuous measure- 
ment of observables that are linear combinations of the 
canonical variables 20]. In the general case of measure- 
ments that are not QND the equivalent classical model 
will have noise-driven dynamical equations as well as 
noise on the measured signal. The noise processes driving 
the dynamics and the measured signal may also be cor- 
related. The case of position measurement of a harmonic 
oscillator shows all of these features 28]. 



APPENDIX B: RICCATI EQUATION SOLUTION 
METHOD 



The matrix Riccati equation is ubiquitous in optimal 
control. Here, following 29], we show how to reduce the 
non-linear problem to a set of linear differential equa- 
tions. Consider the generic Riccati Equation: 



dV{t) 
dt 



C - DV(t) - \{t)A - \{t)B\{t) 



We propose the decomposition: 
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with the hnear dynamics 



~dWit)' 

dt 




-D C 




W{t) 


dUlt) 
dt 




B A 




V{t)_ 



It is straightforward to then show that this hnearized 
solution is equivalent to the Riccati equation 



dV{t) dW{t) 



dt 



dt 



U 



■W{t) 



dv-\t) 



dt 



= (-DW(t) +CU(t))U-^(t) 

-W{t)U-\t){BW{t) + AU(t))U-^(t) 
= C - DV(t) - V{t)A - V{t)BV{t) 

where we have used the identity 

^ = -u-(*)«u-(*) 

dt ^ ^ dt ^ ^ 

Thus the proposed solution works and the problem can 
be solved with a linear set of differential equations. 



APPENDIX C: ROBUST CONTROL IN 
FREQUENCY SPACE 

Here we ap ply traditional frequency-space robust con- 
trol methods [s^, El to the classical version of our sys- 
tem. This analysis is different from the treatment in the 
body of the paper in several respects. First, we assume 
nothing about the noise sources (bandwidth, strength, 
etc.). Also, this approach is meant for steady state sit- 
uations, with the resulting estimator-controller being a 
constant gain transfer function. The performance crite- 
rion we present here is only loosely related to the more 
complete estimation description above. Despite these dif- 
ferences, this analysis gives a very similar design proce- 
dure for the steady state situation. 

We proceed as follows with the control system shown 
in FigureEl where we label h{t) = u{t) -\-b{t) as the total 
field. Consider the usual spin system but ignore noise 
sources and assume we can measure z{t) directly, so that 
z{t) = y(t). For small angles of rotation, the transfer 
function from h{t) to y{t) is an integrator 



dyit) dz{t) 



dt 
sy{s) 

P{s) 



dt 
jJh{s) 

P{s)h{s) 

jj/s 



-iJh(t) 



Now we define the performance criterion. First notice 
that the transfer function from the field to be measured 
h{t) to the total field h{t) is S{s) where 



h{s) 
S{s) 



S{s)b{s) 
1 



l^P{s)C{s) 



(Also notice that this represents the transfer function 
from the reference to the error signal e{s) = S{s)r{s).) 
Because our field estimate will be b{t) = —u{t)^ we de- 
sire h{t) to be significantly suppressed. Thus we would 
like S{s) to be small in magnitude (controller gain \C{s)\ 
large) in the frequency range of interest. However, be- 
cause the gain \C{s)\ must physically decrease to zero at 
high frequencies we must close the feedback loop with 
adequate phase margin to keep the closed-loop system 
stable. This is what makes the design of C{s) non-trivial. 

Proceeding, we now define a function Wi {s) which rep- 
resents the degree of suppression we desire at the fre- 
quency s = juo. So our controller C{s) should satisfy the 
following performance criterion 



\Wi{s)S{s) 



< 1 



Thus the larger Wi{s) becomes, the more precision we 
desire at the frequency s. We choose the following per- 
formance function 



^10 



s/uji 



such that uji is the frequency below which we desire sup- 
pression 1/VFio. 

Because our knowledge of J is imperfect, we need to 
consider all plant transfer functions in the range 



S 



J max\ 



Our goal is now to find a C{s) that can satisfy the 
performance condition for any plant in this family. We 
choose our nominal controller as 



Co (5) 



So if J = J' then the system closes at cjc (i-e., 
\P{iUi)c)CQ{iuoc)\ = 1, whereas in general the system will 
close at cocR = jr • We choose this controller because 
P{s)C{s) should be an integrator (ex 1/s) near the clos- 
ing frequency for optimal phase margin and closed loop 
stability. 

Next we insert this solution into the performance con- 
dition. We make the simplifying assumption coi <^ loc-jt 



Feedback 
Controller 



b(t) 



c 


u(t) + 














Atomic 
Ensemble 



y(t) 



b(t) 



FIG. 8: Spin control system with plant transfer function 
P{s) = iJ/s. r{t) is the reference signal, which is usually 
zero. e(t) is the error signal. u{t) is the controller output. 
b{t) is the external field to be tracked. h{t) = b{t) + u{t) is 
the total field. b{t) is the field estimate. 
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(we will check this later to be self-consistent). Then the 
optimum of the function is obvious and the condition of 
Eq. (|('1|) becomes 



UJiWio < LOCR iOc 



We want this condition to be satisfied for all possible spin 
numbers, so we must have 



10 = m.\n\ucR\ 



■ LOC- 



Jrr 



J' 



(Cl) 



Experimentally, we are forced to roll-off the controller 
at some high frequency that we shall call ujq . Electronics 
can only be so fast. Of course, we never want to close 
above this frequency because the phase margin would 
become too small, so this determines the maximum J 
that the controller can reliably handle 



luq max[cjci?] 



Jyy 



J' 



(C2) 



Combining Eqs. (|(]1|) and (|(]2|) we find our fundamen- 
tal trade-off 



UlWlQ = UJQ 



Jrr 



Jrr 



(C3) 



which is the basic result of this section. Given experimen- 
tal constraints (such as Jmim Jmaxi and 00 q)^ it tells us 
what performance to expect (1/VKio suppression) below 
a chosen frequency uoi. 
From Eq. JHSI, 

we recognize that the controller gain 
at the closing frequency needs to be 



In the final analysis, we do not need to use J' and ujc 
to parametrize the controller, only the trade-off and the 



gain. Also, notice that now we can express min[a;c'i?] 
To check our previous assumption 



J 



^1 < j7 



10 



Jrr 



which is true if ^ 1. 

Finally, the system will never close below the frequency 
min[cjc'i^] so we should increase the gain below a fre- 
quency ujh which we might as well set equal to min[a;c'i^]. 
This improves the performance above and beyond the cri- 
terion above. Of course we will be forced to level off the 
gain at some even lower frequency uj^ because infinite 
DC gain (a real integrator) is unreasonable. So the final 
controller can be expressed as 



C{s) = \C\cz 



1 



^i/(l + s/ljh) 



1 + s/ijQ cc;l(1 + s/cjl) 
with the frequencies obeying the order 

ool < 

ujh = min[cc;ci?] = ojiWiq < 
J 



Jrr 



max[cjci?] 



J rr 



Jrr 



-uJlW^ 
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Notice that the controller now looks like the steady 
state transfer function in Figure [31 derived from the 
steady state of the full dynamic filter. (The notation 
is the same to make this correspondence clear). Here ujq 
was simply stated, whereas there it was a function of A 
that went to infinity as A ^ 00. Here the high gain due to 
uol and uoh was added manually, whereas before it came 
from the design procedure directly. 
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